// STL
#include <algorithm>              // std::max_element
#include <cmath>                  // std::pow()
#include <functional>             // std::functional
#include <vector>                 // std::vector

// stats++
#include "statsxx/statistics.hpp" // statistics::mean(), ::stddev_p()

// this
#include "statsxx/optimization/EA/EAParam.hpp"    // EAParam
#include "statsxx/optimization/EA/Individual.hpp" // Individual
#include "statsxx/optimization/EA/Population.hpp" // Population
#include "statsxx/optimization/EA/utility.hpp"    // calc_genome_distance()


//
// DESC: Calculates the shared fitness of an individual with a population.
//
// NOTES:
//     ! The fitness sharing function was taken from:
//          S. W. Mahfoud, PhD thesis
//
inline double calc_shared_fitness(
                                  const Individual &indiv,
                                  const Population &population,
                                  // -----
                                  const EAParam &param
                                  )
{
    // AT LEAST 3 POINTS ARE NEEDED TO CALCULATE A STANDARD DEVIATION IN ORDER TO CALCULATE A SHARING FUNCTION BELOW ...
    if( population.individuals.size() < 3 )
    {
        return indiv.fitness;
    }


    // CALCULATE HOW WELL GENOMES ARE SEPARATED IN THE POPULATION, AND DETERMINE A MIN. SEPARATION TO CATEGORIZE TWO DISTANCES AS THE SAME NICHE ...
    std::vector<double> sep;
    for( decltype(population.individuals.size()) i = 0; i < (population.individuals.size()-1); ++i )
    {
        for( decltype(population.individuals.size()) j = (i+1); j < population.individuals.size(); ++j )
        {
            sep.push_back( calc_genome_distance(population.individuals[i].genome, population.individuals[j].genome) );
        }
    }

    // NOTE: We (obviously) have the total population available; so calculate the population standard deviation.
    double dmin = statistics::mean(sep) - 2.0*statistics::stddev_p(sep);


    // NOW CALCULATE SHARED FITNESS ...
    double sh = 1.;
    for( auto &other : population.individuals )
    {
        double d = calc_genome_distance(indiv.genome, other.genome);

        if( d < dmin )
        {
            sh += (1. - std::pow( (d/dmin), param.fitness_share_alpha ));
        }
    }


    return (indiv.fitness/sh);
}
